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Abstract: The development of benchmark examples for quasi-static delamination 
propagation and cyclic delamination onset and growth prediction is presented and 
demonstrated for Abaqus/Standard. The example is based on a finite element model of a 
Double-Cantilever Beam specimen. The example is independent of the analysis software 
used and allows the assessment of the automated delamination propagation, onset and 
growth prediction capabilities in commercial finite element codes based on the virtual 
crack closure technique (VCCT). First, a quasi-static benchmark example was created 
for the specimen. Second, based on the static results, benchmark examples for cyclic 
delamination growth were created. Third, the load-displacement relationship from a 
propagation analysis and the benchmark results were compared, and good agreement 
could be achieved by selecting the appropriate input parameters. Fourth, starting from 
an initially straight front, the delamination was allowed to grow under cyclic loading. 
The number of cycles to delamination onset and the number of cycles during 
delamination growth for each growth increment were obtained from the automated 
analysis and compared to the benchmark examples. Again, good agreement between the 
results obtained from the growth analysis and the benchmark results could be achieved 
by selecting the appropriate input parameters. The benchmarking procedure proved 
valuable by highlighting the issues associated with choosing the input parameters of the 
particular implementation. Selecting the appropriate input parameters, however, was not 
straightforward and often required an iterative procedure. Overall the results are 
encouraging, but further assessment for mixed-mode delamination is required. 
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1. introduction 

The virtual crack closure technique (VCCT) (Rybicki and Kanninen, 1977; Krueger, 2004) is 
widely used for computing energy release rates based on results from finite element analyses and 
to supply the mode separation required when using a mixed-mode fracture criterion to predict 
delamination onset or growth. As new methods for analyzing composite delamination are 
incorporated in finite element codes, the need for benchmarking becomes important since each 
code requires specific input parameters unique to its implementation. Once the parameters have 


2012 SIMULIA Customer Conference 


1 



been identified, they may then be used with confidence to model delamination growth for more 
complex configurations. 

An approach for assessing the mode I, mode II and mixed-mode I/II, delamination propagation 
capabilities in commercial finite element codes was recently presented and demonstrated for 
Abaqus/Standard (Krueger, 2008, 2010, 2011, 2012) as well as MD Nastran™ and Marc™ 2 (Orifici 
and Krueger, 2010). The current paper summarizes the development of benchmark cases based on 
the Double Cantilever Beam (DCB) specimen, and the creation of examples is shown step-by-step. 
The creation step is independent of the software used. After creating the benchmark, the approach 
is demonstrated for the automated analysis in the commercial finite element code 
Abaqus/Standard. Starting from an initially straight front, the delamination is allowed to grow 
based on the algorithms implemented into the commercial finite element software. Input control 
parameters are varied to study the effect on the computed delamination propagation and growth. 
Results obtained from the automated analysis should closely match the results created manually to 
obtain the benchmark example. 

2. Specimen and Model Description 


For the current numerical investigation, the Double Cantilever Beam (DCB), as shown in Figure 1, 
was chosen. The DCB specimen is used to determine the mode I interlaminar fracture toughness, 


Gic- 



dimensions 

B 25.0 mm 

2h 3.0 mm 

2L 150.0 mm 

a 0 30.5 mm 

layup: [0] 24 

static loading 

6/2 2.0 mm 

fatigue loading 

6 max /2 0.67 mm 
S min /2 0.067 mm 
R 0.1 

f 10.0 Hz 


Figure 1. Double Cantilever Beam (DCB) Specimen. 


Typical two- and three-dimensional finite element models of the DCB specimen are shown in 
Figures 2 and 3. Along the length, all models were divided into different sections with different 
mesh refinement. The specimen was modeled with six elements through the specimen thickness. 
The plane of delamination was modeled as a discrete discontinuity in the center of the specimen. 


2 MD Nastran™ and Marc™ are manufactured by MSC. Software Corp., Santa Ana, CA, USA. 
NASTRAN®’ is a registered trademark of NASA. 
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To create the discrete discontinuity, each model was created from separate meshes for the upper 
and lower part of the specimens with identical nodal point coordinates in the plane of 
delamination. Two surfaces (top and bottom surface) were created on the meshes as shown in 
Figure 3. Additionally, a node set was created to identify the intact (bonded nodes) region. The 
initial delamination front is highlighted in red. The modeling is discussed in detail in related 
reports (Krueger, 2008, 2010). 



Figure 2. Deformed model of a DCB specimen (plane strain-CPE4). 



bonded 

nodes 


top and bottom 
surface 


initial straight delamination front 


Figure 3. Full three-dimensional (solid C3D8I) finite element model of a DCB 

specimen. 
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3. Creating a Benchmark Solution for Quasi-Static Loading 


First, models simulating specimens with seven different delamination lengths (30 mm £ a 0 £ 
40mm) were analyzed. For each delamination length modeled, the reaction loads P at the location 
of the applied displacement were calculated and plotted versus the applied opening displacement 
8/2 as shown in Figure 4. 



Figure 4. Calculated critical behavior and resulting benchmark case obtained from 
analyses using three-dimensional models (C3D8I). 

The mode I strain energy release rate was calculated along the delamination front across the width 
of the specimen. The critical load, P crit , when the failure index in the center of the specimen 
( y/B=0 ) reaches unity ( Gj/G c =l ), can be calculated based on the relationship between load P and 
the energy release rate. 


P 2 dC P 
2 dA 


( 1 ) 


In Equation (1), Cp is the compliance of the specimen and dA is the increase in surface area 
corresponding to an incremental increase in load or displacement at fracture. The critical load P cri , 
and critical displacement 8 crit /2 were calculated for each delamination length modeled. 
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and the results were included in the load/displacement plots as shown in Figure 4 (solid red 
circles). The results indicate that, with increasing delamination length, less load is required to 
extend the delamination. This means that the DCB specimen exhibits unstable delamination 
propagation under load control. Therefore, prescribed opening displacements d/2 were applied in 
the analysis instead of nodal point loads P to avoid problems with numerical stability of the 
analysis. The computed critical load/displacement results were used as a benchmark. The 
benchmark result may also be visualized by plotting the prescribed opening displacements d/2 at 
delamination growth onset versus the delamination length, a, as illustrated in Figure 5. For the 
delamination propagation, therefore, the load/displacement results obtained from the model of a 
DCB specimen with an initially straight delamination of a=30 mm length should closely match the 
critical load/displacement path (solid red line) in Figure 4 or the prescribed opening 
displacements/delamination-length relationship Figure 5. 



Figure 5. Critical delamination length-displacement behavior for DCB specimen 
obtained from analyses using three-dimensional models (C3D8I). 
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4. Results from Automated Propagation Analyses 


A set of example results are shown in Figure 6 where the computed resultant force (load P) at the 
tip of the DCB specimen is plotted versus the applied crack tip opening (<5/2) for different input 
parameters. 
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applied opening displacement <5/2, mm 

Figure 6. Computed load-displacement behavior obtained from automated 
analysis using three-dimensional models (C3D8I). 

To overcome the convergence problems, the methods implemented in Abaqus/Standard were used 
individually to study the effects. First, global stabilization was added to the analysis. For a 
stabilization factor of 2x1 0" 5 , the stiffness changed to almost infinity once the critical load was 
reached causing the load to increase sharply (plotted in red). The load increased until a point was 
reached where the delamination propagation started and the load gradually decreased following a 
saw tooth curve with local rising and declining segments. The gradual load decrease followed the 
same trend as the benchmark curve (in grey) but is shifted toward higher loads. For a stabilization 
factor of 2xl0" 6 (in green), the same saw tooth pattern was observed but the average curve was in 
good agreement with the benchmark result. Second, contact stabilization was added to the 
analysis. For all combinations of stabilization factors and release tolerances, a saw tooth pattern 
was observed, where the peak values were in good agreement with the benchmark result (example 
plotted in blue). Third, viscous regularization was added to the analysis to overcome 
convergence problems. Convergence could not be achieved over a wide range of viscosity 
coefficients when a default release tolerance value of 0.2 was used. Subsequently, the release 
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tolerance value was increased. For all combinations of the viscosity coefficient and release 
tolerance, where convergence was achieved, a saw tooth pattern was obtained, where the peak 
values were in good agreement with the benchmark result (in black). 

An alternative way to plot the benchmark is shown in Figure 7 where the prescribed opening 
displacements 812 at delamination growth onset is plotted versus the delamination length, a, This 
way of presenting the results is shown, since it may be of advantage for large structures where 
local delamination propagation may have little influence on the global stiffness of the structure 
and may therefore not be visible in a global load/displacement plot. However, extracting the 
delamination length a from the finite element results required more manual, time consuming post- 
processing of the results compared to the relatively simple and readily available output of nodal 
displacements and forces. The results plotted in Figure 7 are selected examples that were 
discussed above. The conclusions that can be drawn from this plot are identical to those discussed 
above. 



Figure 7. Computed delamination length-displacement behavior obtained from 
automated analysis using three-dimensional models (C3D8I). 

The examples shown highlight the importance of benchmarking to identify critical analysis input 
parameters. In summary, good agreement between the load-displacement relationship obtained 
from the propagation analysis results and the benchmark results could be achieved by selecting the 
appropriate input parameters. However, selecting the appropriate input parameters such as release 
tolerance, global or contact stabilization and viscous regularization, was not straightforward and 
often required an iterative procedure. The default setting for global stabilization yielded 
unsatisfactory results although the analysis converged. The detailed results for Abaqus/Standard 
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are discussed in detail in a related report (Krueger, 2008). The results of a similar study using the 
same approach for the finite element codes Marc and MSC.Nastran were published recently 
(Orifici and Krueger, 2010). 

5. Creating a Benchmark Solution for Cyclic Loading 

For the cyclic loading of the specimen, guidance was taken from a draft standard designed to 
determine mode I fatigue delamination propagation. In the draft document, it is recommended to 
start the test at a maximum displacement, <5 max , which causes the energy release rate at the front, 
Gimax, to reach initially about 80% of G Ic . For the current study, a critical energy release rate 
G/ c =0.17 kJ/m 2 was used and the critical values P crit and 5 crit (grey dashed lines) were obtained 
from the benchmark for static delamination propagation shown in the load-displacement plot in 
Figure 8. 
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Figure 8. Computed load-displacement behavior obtained from automated 
analysis using three-dimensional models (C3D8I). 

The maximum load, P max , and maximum displacement, d max /2, were calculated using the known 
quadratic relationship between energy release rate and applied load or displacement given in 
Equation (2). The calculated maximum load, P max , and calculated maximum displacement, d max /2, 
are shown in Figure 8 (dashed red line) in relationship to the static benchmark case (solid grey 
circles and dashed grey line) mentioned above. During constant amplitude cyclic loading of a 
DCB specimen under displacement control, the applied maximum displacement, d max /2=0.61 mm, 
is kept constant while the load drops as the delamination length increases (solid red circles and 
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solid red line). For each of die seven finite element models representing seven delamination 
lengths (30.0 mm< a fl <40.0 mm) as shown in Figure 4, the energy release rate corresponding to an 
applied maximum displacement d max l2=0.61 mm was calculated using Equation (2). The energy 
release rate decreases with increasing delamination length, a, as shown in Figure 9 (solid red 
circles and solid red line). Delamination growth was assumed to stop once the calculated energy 
release rate drops below the cutoff value, G,h, (green solid horizontal line). The static benchmark 
case (solid grey circles and dashed grey line in Figure 9), where the delamination propagates at 
constant Gj c (solid blue line in Figure 9) was included for comparison. 



Figure 9. Energy release rate - delamination length behavior for a DCB specimen 
obtained from two-dimensional models (CPE4). 


The number of cycles to delamination onset, No, can be obtained from the delamination onset 
curve which is a power law fit of experimental data obtained from a DCB test using the respective 
standard for delamination growth onset as shown in Figure 10. Solving for No yields 

l i 


G = m 0 ■ N D ‘ 


N d = 


Nd - c r G Cl , c 2 - 


( 3 ) 


\ "'o 


At the beginning of the test, the specimen is loaded initially so that the energy release rate at the 
front, G Imax , reaches about 80% of G Ic corresponding to G /max =0 . 1 362 kj/m 2 . From the 
delamination onset curve, the number of cycles to delamination onset is determined as, Nd= 150. 
Details are discussed in a related report (Krueger, 2010). 
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Figure 10. Experimentally determined delamination growth onset for DCB 

specimen. 


The number of cycles during stable delamination growth, Nc, can be obtained from the fatigue 
delamination propagation relationship (Paris Law) as shown in Figure 11. The delamination 
growth rate (solid purple line) can be expressed as a power law function 


da_ 

dN 


or 


A a 
AN 


G n 
n 


(4) 


where da/dN is the increase in delamination length per cycle and G max is the maximum energy 
release rate at the front at peak loading. The factor c and exponent n are obtained by fitting the 
curve to the experimental data obtained from DCB tests. A cutoff value, G t h, was chosen below 
which delamination growth was assumed to stop. Details are discussed in a related report 
(Krueger, 2010). Note, that this benchmarking exercise ignores branching, or fiber bridging. 
Hence, the Paris Law was not normalized with the static R-curve as recently suggested (Krueger, 
2010 ). 

For the current study, increments of Aa= 0.1 mm were chosen. Starting at the initial delamination 
length a 0 =30.5 mm, the energy release rates G Unax were obtained for each increment, i, from the 
curve fit (solid red circles and solid red line) plotted in Figure 9. These energy release rate values 
were then used to obtain the increase in delamination length per cycle or growth rate Aa/AN from 
the Paris Law. The number of cycles during stable delamination growth, N c „ was calculated by 
summing the increments AN-, 
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( 5 ) 


N c =^ AN i = ^-G;" max - A a => a = a 0 + 

i-1 i-\ C i-l 

where k is the number of increments. The resulting delamination length, a, was calculated by 
adding the incremental lengths Aa to the initial length a 0 . 


^ Aa = a Q + k ■ Aa 


da/dN, 

mm/cycle 



Figure 11. Experimentally determined delamination growth rate (Paris Law). 

For the combined case of delamination onset and growth, the total life, Nr, may be expressed as 
N t = N d +N c where, N», is the number of cycles to delamination onset and N G , is the number of 
cycles during delamination growth. For this combined case, the delamination length, a, is plotted 
in Figure 12 for an increasing number of load cycles Nt- For the first N n cycles, the delamination 
length remains constant (horizontal red line), followed by a growth section where - over N G cycles 
- the delamination length increases following the Paris Law (crosses and solid red line). Once a 
delamination length is reached where the energy release rate drops below the assumed cutoff 
value, G,h, (as shown in Figure 7), the delamination growth no longer follows the Paris Law 
(dashed grey line) and stops (horizontal solid red line). A delamination length prediction analysis 
that accounts for delamination fatigue onset as well as stable growth should yield results that 
closely resemble the plot in Figure 12. The curve fit (solid red line) can therefore be used as a 
benchmark. 
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Figure 12. Delamination onset and growth behavior for a DCB specimen obtained 
from two-dimensional models (CPE4). 


6. Results from Automated Growth Analyses 


A set of example results are shown in Figure 13 where the delamination length, a, is plotted versus 
the number of cycles, N, for different input parameters and models. For all results shown, the 
analysis stopped when a 10,000,000 cycles limit - used as input to terminate the analysis - was 
reached. 

First, the effect of the initial time increment used in the analysis settings was studied as shown in 
Figure 13. The initial time increment was varied between / fl =0.01 (one tenth of a single loading 
cycle, f,=0.1s, open blue circles and solid blue line) and z' 0 =0.0001 (one thousandth of a single 
loading cycle, open green diamonds and solid green line). For larger initial time increments, the 
onset of delamination shifted towards a lower number of cycles. Reducing the initial time 
increments, however, significantly increased the computation time. Based on the results, it was 
therefore decided to use an initial time increment of z'o=0.001 (open red squares and solid red line) 
for the remainder of the study to save computation time. This step is justified by the fact that the 
results obtained for / fl =0.001 were almost identical to the values obtained from the analysis where 
a smaller initial time increment was used (/ 0 =0.0001). Second, the release tolerance was varied. In 
the current study, varying the input between 0.2 (the default value) and 0.01 did not have any 
effect on the computed onset and growth behavior during the analysis. It is assumed that the 
release tolerance only affects static delamination propagation and not growth under cyclic loading 
studied here. Stabilization (to overcome the convergence problems) which was required for the 
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automated delamination propagation under static loading, was not required for automated growth 
analysis under cyclic loading. However input parameters to define the cyclic loading were varied 
and showed different affects on the computed results. Also, solution controls had to be 
manipulated to achieve reasonable computation times. The results from these studies are discussed 
in detail in a related report (Krueger, 2010). 


delamination 
length 
a, mm 



Figure 13. Computed delamination onset and growth obtained from two- 
dimensional models (CPE4) for different initial time increments i 0 . 

7. Conclusions 

The development of two benchmark examples for static delamination propagation and cyclic 
delamination growth prediction were presented and demonstrated for the commercial finite 
element code Abaqus/Standard. The example was based on a finite element model of a Double 
Cantilever Beam (DCB) specimen, which is independent of the analysis software used and allows 
the assessment of the delamination propagation and growth prediction capabilities in commercial 
finite element codes. The development of the benchmark examples was presented step by step. 
Starting from an initially straight front, the delamination was then allowed to propagate or grow 
under cyclic loading using automated analyses. 

With respect to benchmarking, the results showed the following: 

• Benchmark examples for the assessment of automated delamination propagation and 
growth capabilities can easily be developed from a set of static analyses of models with 
different initial delamination lengths. 
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• The development of benchmark examples is independent of the software used and 
independent of experimental results. 

• The benchmarking process can be used to identify the issues associated with the input of 
a particular finite element code or implementation. 

• Benchmarking helps identify relevant input parameters so that they may be used with 
confidence when modeling more complex configurations. 

In general, good agreement between the results obtained from the propagation and growth analysis 
and the benchmark results could be achieved by selecting the appropriate input parameters. 
Overall, the results are promising. In a real case scenario, however, where the results are unknown, 
obtaining the right solution will remain challenging. Further studies are required which should 
include the assessment of the propagation capabilities in more complex mixed-mode specimens 
and on a structural level. 
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